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Abstract. We show that multi-orbital and density-induced tunneling have a significant 
impact on the phase diagram of bosonic atoms in optical lattices. Off-site interactions lead to 
density-induced hopping, the so-called bond-charge interactions, which can be identified with 
an effective tunneling potential and can reach the same order of magnitude as conventional 
tunneling. In addition, interaction-induced higher-band processes also give rise to strongly 
modified tunneling, on-site and bond-charge interactions. We derive an extended occupation- 
dependent Hubbard model with multi-orbitally renormalized processes and compute the 
corresponding phase diagram. It substantially deviates from the single-band Bose-Hubbard 
model and predicts strong changes of the superfluid to Mott-insulator transition. In general, 
the presented beyond-Hubbard physics plays an essential role in bosonic lattice systems and 
has an observable influence on experiments with tunable interactions. 
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1. Introduction 

Hubbard models are extremely successful in describing a variety of systems ranging from 
electrons in solids to bosonic quantum gases in optical lattices. While in solids all kinds of 
complex lattice systems are realized, optical lattices offer simple geometries in combination 
with precisely controllable experimental conditions. The striking idea that ultracold bosonic 
atoms in optical lattices constitute an excellent tool for studying the superfluid to Mott- 
insulator transition [1] has initiated an intensive investigation of bosonic lattice models. It 
has caused a large number of theoretical studies on the Bose-Hubbard model ranging from 
mean-field treatments [2-4] and quantum Monte Carlo [5] to DMRG [6] approaches. In 
addition, the influence of disorder [2, 7-11] and the effect of confining potentials have been 
discussed [12-15]. However, the first-order corrections to the Bose-Hubbard model itself, 
namely higher orbital tunneling (Fig. la) and bond-charge interactions (Fig. lb), have so 
far been mostly neglected in these approaches. Bond-charge interactions known from solids 
[16-19] constitute a density-induced tunneling process caused by the scattering of particles 
on neighboring sites. This extension to the Hubbard model was discussed in the context of 
superconductivity [16] and ferromagnetism [19]. While for fermionic systems these effects 
are partly suppressed due to the Pauli principle, they have large impact on bosonic systems. 
The unique experimental access in optical lattices allows for a detailed investigation of this 
interaction-assisted tunneling. This density-dependent hopping was recently also discussed 
for boson-fermion mixtures [20, 21], where the Bose-Fermi-Hubbard model does not fully 
cover effects arising from the interspecies interaction [20-24]. In this context, the role of 
higher orbitals for the tunneling processes [20-22] was investigated. For bosonic systems, 
modifications of tunneling and on-site interaction due to higher bands were studied by 
variational mean-field methods [25-28] and by numerical exact methods restricted either to 
double- or triple- well systems [29-32] or to the on-site interaction [33-36]. However, in 
the correlated regime present at the vicinity of the Mott-insulator transition the applicability 
of mean-field methods to describe higher-band processes is doubtful. Thus, the challenging 
problem is to find a comprehensive description including both bond-charge hopping [16- 
20, 37, 38] and higher band processes [25-28, 30, 31, 33-36, 39, 40] that is valid also for 
strongly correlated systems. The recent experimental progress in accessing the occupation- 
dependent on-site interactions [35, 41-43] and the tunneling matrix element [24] requests for 
an accurate theoretical method to calculate these parameters. Furthermore, the superfluid to 
Mott-insulator transition in experiments with tunable interactions [23, 42] is directly affected 
by bond-charge tunneling as well as multi-orbital extensions to conventional and bond-charge 
tunneling. 

Here, we show that interaction-induced processes cause two substantial modifications of 
the Bose-Hubbard model for bosonic atoms in three-dimensional optical lattices. First, so- 
called bond-charge or density-induced tunneling combining interaction and hopping (Fig. lb) 
is a notable contribution to tunneling. Second, higher bands are occupied due to the interaction 
of particles on a lattice site. This increases the tunneling significantly as the tunneling 
in higher bands (Fig. la) is strongly enhanced. We present a general scheme to calculate 
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Figure 1: Hopping processes beyond the standard single-band Hubbard tunneling: Example 
processes for (a) multi-orbital tunneling, (b) single-orbital bond-charge assisted hopping and 
(c) multi-orbital bond-charge hopping, (d) Initial and final states for one particle tunneling 
from the left to the right site. The interaction causes an admixture of higher orbital states 
requiring a renormalization of the bare lowest-band tunneling processes. 



renormalized tunneling processes by effectively incorporating higher-band contributions. On 
this basis, we derive an effective occupation-dependent model where the renormalization 
solely depends on the solution of the on-site problem. The well-known phase diagram for 
bosons in lattices describing the superfluid (SF) to Mott insulator (MI) transition is drastically 
altered. The presented results have direct relevance for optical lattice experiments, which can 
therefore serve as an ideal testing ground for beyond Hubbard physics. 

After presenting the resulting phase diagram in section 2, we discuss the extensions of 
the Hubbard model individually starting with the effect of bond-charge interactions in section 
3. The multi-orbital renormalization of the operators for conventional tunneling, bond-charge 
tunneling and on-site interaction is elaborated in sections 4-7. The basic idea is to separate the 
strongly-correlated on-site problem from the inter-site dynamics. First, we solve the on-site 
problem numerically exactly using the Wannier states of the lattice. Subsequently, we derive 
renormalized matrix elements for the multi-orbitally dressed operators. Finally, we calculate 
the phase diagram by treating the many-site problem with mean-field theory in section 8. Note 
that the multi-orbital renormalization can be used for various methods to treat the many-site 
problem and that the results are expected to be similarly affected. 



2. Phase diagram and model Hamiltonian 

We now first explain as one central result the phase diagram for bosonic atoms (Fig. 2) 
which is derived and further discussed afterwards. The black lines depict the Mott lobes 
predicted by the single-band Bose-Hubbard model applying mean-field theory [4]. Our 
calculations (colored lines) take into account the corrections to the tunneling as mentioned 
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Figure 2: Phase diagrams showing the superfluid to Mott- insulator transition for interaction 
strengths a s normalized to the lattice spacing a (a) a s /a = 0.042 and (b) a s /a = 0.014. The 
latter corresponds to the background scattering length a s = 100 a of 87 Rb at a lattice spacing 
of a = 377 nm [23]. The phase boundaries are plotted for the Bose-Hubbard model (black) 
and the full Hamiltonian Eq. (1) using mean-field (blue) and Gutz wilier (dashed red). The 
green line corresponds to the occupation-independent Hamiltonian Eq. (2). 



above as well as multi-orbital modifications of the on-site interaction. The phase diagrams 
are obtained for bosonic atoms of arbitrary species in cubic sinusoidal optical lattices with 
V(r) = Vo J2i cos 2 (7rr^/a), where i = {x, y, z], a is the lattice spacing and Vo the lattice 
depth. We predict that the phase area of the Mott insulator is considerably reduced even 
for moderate interaction strengths and substantially deformed for stronger interactions. This 
corresponds to a significant shift of the critical lattice depth of the SF-MI transition (Fig. 3) 
due to an effectively increased tunneling and reduced on-site interaction. The discrepancy can 
reach more than 5 E R forn = 3 particles per lattice site. We expect the shift to be observable 
in experiments with a filling n > 2 and strong or tunable interactions. 

The mentioned extension to the Hubbard model leading to the phase diagram in Fig. 2 
(red and blue lines) can be included in an effective Hamiltonian 

Hm = -£ 6& Jgk + \J2 U^Hiifk ~ I)- (1) 

with multi-orbitally renormalized creation and annihilation operators b\ and 6 i5 respectively, 
and the operator hi = bjft- corresponding to the number of particles on site i. This Hamiltonian 
has the same structure as the Bose-Hubbard Hamiltonian [1, 2] but with multi-orbitally 
renormalized tunneling matrix elements and on-site interactions U~£° that explicitly 

depend on the site occupation n*. In the following, we calculate these occupation-dependent 
parameters with exact diagonalization and derive the phase boundaries using a mean-field 
(MF) and a Gutzwiller (GW) approach. Furthermore, we present a simplified Hamiltonian 
with occupation-independent parameters (green lines in Fig. 2) which is in good agreement 
with the full model above. This Hamiltonian is restricted to conventional and bond-charge 
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Figure 3: Critical lattice depths of the superfluid to Mott- insulator transition for the Hubbard 
model (dashed), the full multi-orbital Hamiltonian Eq. (1) and the occupation-independent 
lowest-band Hamiltonian Eq. (2). The results are plotted as a function of the scattering length 
a s /a in units of the recoil energy E R = h 2 /8ma 2 . 




tunneling in the lowest-band 

^occ.-indep. = - [«/ b\bj + Jbc + + y ^ ^(fii - 1) (2) 

with single-band operators SJ and 6$, the Hubbard tunneling J, lowest-band bond-charge 
interaction J BC and renormalized on-site interaction U 2 for two particles, as is elaborated 
in detail below. In general, the results for Hamiltonian (2) show that the lowest-band bond- 
charge tunneling is an important contribution and its influence increases with the number 
of particles per site. Note that the mean-field and Gutzwiller approach fully coincide 
for occupation-independent models (green and black lobes), whereas there is a minimal 
discrepancy for the Hamiltonian (1) with occupation-dependent tunneling. In the mean-field 
calculation, the decoupling of the lattice sites requires the use of J^ n as an approximation 
for JnXi.rn and J nL±i ( see Appendix F and G). 



3. Bond-charge interactions 

The two-particle interaction of bosonic particles in the lowest band of the lattice is given by 
Ant = \ Ysijki U ijki b\b]hbi with interaction integrals 

Uijki =9 Jd 3 r w*(r)wj(r)w k (r)wi(r) (3) 

and Wannier functions i^(r) = w^(r — r^) of the lowest band, which are maximally 
localized at site i. Here, a ^-interaction potential with repulsive interactions g = ^a s > 
is assumed. One of the key ingredients of the Hubbard model [1, 2, 44] is to restrict the 
two-particle interaction to the dominating term, namely, the on-site interaction U = Una. 
However, some of the neglected off-site processes correspond to scattering accompanied by 
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hopping to a neighboring site (e.g., b\bjbjbj) known as bond-charge tunneling (Fig. lc). It 
is essential that this process constitutes a hopping process which physically modifies the 
overall tunneling in the system. Thus, although this density-dependent hopping contribution 
is small in comparison with the on-site interaction U, it can only be neglected if also small in 
comparison with the conventional tunneling J. Considering two neighboring sites i and j 9 the 
density-induced tunneling operator takes the form 

Jbc = -Jbc b\{fii + hj)bj with J BC = -U m (4) 

and J B c > 0. Due to the structure of this operator, we can combine it with the conventional 
tunneling J to an effective hopping 

=-[J+ Jbc(^ + n 3 ~ 1)] b\bj. (5) 

As a consequence, the total tunneling increase is of the order of several ten percent for standard 
87 Rb-conditions (see Fig. 4, dashed line) and even more for stronger interactions and deeper 
lattices. Note that density-density off-site interactions (£%j) [6] and correlated pair tunneling 
(Uujj) are considerably weaker than the bond-charge interactions (see Appendix C). 

The nature of bond-charge interactions can be intuitively illustrated as an effective 
potential. In Eq. (5), the term hi + h j — 1 corresponds to the density Pbc(f) = ^i\^i{^)\ 2 + 
(rij — l)\wj(r)\ 2 on sites i and j excluding the hopping particle. By inserting the explicit 
expressions for the integral J B c and the tunneling J, we can write the effective hopping 
operator (5) as 

J e B ff c = J d 3 r w* ^ + V(v) + ^Bc(r)) Wj b%. (6) 

We can now identify the expression V(r) + gp BC (r) as an effective tunneling potential which 
is illustrated in Fig. 4c. Since the density of the repulsively interacting particles is localized 
at the centers of the lattice sites, the effective potential corresponds to a shallower lattice, 
i.e., increased tunneling. Note that the effective potential can be used to determine the total 
tunneling via band structure calculations (see Appendix E). 

4. Multi-orbital tunneling 

While the extension to the Hubbard model above is still within the lowest-band 
approximation, the following discussion covers the influence of higher orbitals. The 
calculation of the effective tunneling can be performed by taking two adjacent sites of the 
lattices into account (Fig. Id). Writing initial and final states as a product of wave functions 
of the two sites allows to express the renormalized tunneling in terms of single-site expectation 
values. Consequently, the renormalized parameters depend only on the solution of the 
single-site problem. The presented scheme allows to derive the effective Hamiltonian (1) 
independently of the applied method for treating the on-site many-particle problem. Here, the 
orbital occupations are computed by exact diagonalization of n particles on a single site using 
a finite many-particle basis with a high-energy cut-off. The short-range interaction between 
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Figure 4: Contributions to the effective tunneling J%* n . with rii = n j = 3 by multi-orbital 
tunneling J^° n ., bond-charge interactions J B c and multi-orbital bond-charge interactions 
Jn^nf°> where the latter two scale with the pref actor v = + rij — 1. The deviations 
from the bare Hubbard tunneling J are plotted as a function of (a) the lattice depth V 
and (b) the interaction strength a s . (c) The single-band bond-charge interaction gives rise 
to an effective tunneling potential (see text), (d) Qualitatively, the multi-orbital tunneling 
scales with the population of higher orbitals (shading) where the tunneling J a (arrows) is 
substantially enhanced in higher bands. 



the atoms is modeled as a box-shaped interaction potential (see Appendix A for details). The 
single-particle orbitals are represented by Wannier functions (r) with a three-dimensional 
band index a = {oi X) a y ,a z ) . 

In the following, we consider multi-orbital hopping processes, where one atom tunnels 
from the left site (L) with initially n L particles to the right site (R) with initially n R 
particles (Fig. Id). As a direct consequence, the total tunneling J nL ,n R becomes intrinsically 
occupation-dependent. The multi-orbital hopping operator J MO = J2 a J<* b^b^ sums over 
all possible multi-orbital processes 

Ja = -J d 3 r w£> + V(r)j w£\ (V) 

where tunneling takes place only between equal orbitals due to orthogonality relations. It is 
possible to reduce this complex multi-orbital operator to an effective hopping operator 

j£° = -%hJ%k. (8) 

The occupation-dependent tunneling matrix elements are given by J^ nR oc (^ F | J MO 
where = ^(n L , n R ) denotes the initial and \I> F = \I/(n L — 1, n R + 1) the final state (Fig. Id). 
This expression can be separated into matrix elements for the respective lattice sites 

Cln = 7 = ^ E J « + h(a)] l*M> ^ - ^ l*M> ' (9) 

V n L(n R + 1) ^ 

where ty(n) is the single-site wave function. This can be evaluated using the coefficients 
obtained by exact diagonalization as elaborated in Appendix B. 
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The tunneling in higher bands J a is drastically enhanced compared with the lowest band 
tunneling J . Therefore, the contributions of higher-orbital processes are significant and grow 
exponentially with the lattice depth, although the population of higher bands is usually smaller 
than 1%. As an example, the matrix element J^ 3 ° for three bosons per site is plotted in 
Fig. 4 (green lines). In general, the effective tunneling is increased except for the elements 
J^q = Jf^n-i where the tunneling is slightly reduced. 

In an intuitive picture, the many-particle calculation can be reduced to the occupation of 
single-particle orbitals (h^). The fraction of particles occupying higher orbitals a tunnels via 
the respective tunneling matrix elements J a , as indicated in Fig. 4(d). In this approximation, 
the effective tunneling can be estimated by the weighted sum (see Appendix B) 

« R - - / = = E W^k^U < 10 > 

This simplified approach is in qualitative agreement with the fully correlated calculation 
above. 

5. Multi-orbital bond-charge interaction 

It proves necessary to apply the same multi-orbital treatment of the normal tunneling 
developed in the last section to the bond-charge assisted hopping. By analogy, we can define 
an effective multi-orbital bond-charge hopping 

j e ^ M ° =-(,!(», +fi>,<sr° (id 

with an occupation-dependent parameter. It depends on the single-site matrix elements 
- 1)| |*(n)) and (*(n - 1)| ti^b^U^ |*(n)) as discussed in Appendix D. The 
results shown in Fig. 4 for the multi-orbital bond-charge interaction (solid red line) differ 
substantially from the single-orbital case (dashed line) for intermediate and deep lattices. 
Processes involving higher bands start to dominate for deep lattices, rendering the lowest- 
band calculation invalid, and can even turn the sign of the bond-charge contributions negative. 
Processes of the type J n ' are only weakly influenced by multi-orbital corrections. 

6. Effective tunneling 

The total tunneling consists of normal and bond-charge assisted tunneling, both effectively 
including higher orbital processes. The total tunneling 

■C = JnX + (ni+nj-l)^ d2) 

and the individual contributions are shown in Fig. 4a and b. The occupation-dependent 
enhancement of the tunneling is depicted in Fig. 5. The deviation from the Hubbard tunneling 
J can easily reach 30% for three particles per site at moderate a s and Vo. The data in this 
figure is calculated using a box-shaped interaction potential (see Appendix A) with a width 
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Figure 5: Occupation-dependent total tunneling J^ n . and on-site interactions as a 

function of (a) the lattice depth Vo and (b) the interaction strength a s . The empirical fit U 2 is 
plotted as a dashed black line. Note that J nj ,m and J ni +i^ nj -i describe time-reversed hopping 
processes and thus are equivalent. For an explanation of the error bars see text. 

W = 5 nm and a lattice constant a = 377 nm. It is important to note that the results are 
almost independent of W, providing that W is much smaller than the lattice constant. The 
error bars in Fig. 5 show the results for W = 25 nm and W — >> 0. 

While the normal tunneling is increased by the multi-orbital renormalization, the density- 
induced tunneling is reduced. Coincidentally, both multi-orbital corrections compensate each 
other and the overall deviations can be approximated surprisingly well by the single-band 
bond-charge tunneling for moderate particle numbers (dashed line in Fig. 4). This justifies the 
applicability of the model Hamiltonian (2), where tunneling is restricted to the lowest band, 
i- e -> Jn% ~J+{n l + n 3 - 1) J BC . 

7. Occupation-number-dependent on-site interaction 

In a similar fashion as the inclusion of higher orbitals alters the tunneling, also the on-site 
interaction U is modified. The orbital degree of freedom decreases the on-site interaction [34- 
36] as the particles tend to avoid each other. The results for U^° 9 representing the eigenvalue 
of the single-site exact diagonalization [35], are depicted for occupation numbers n = 2-4 
in Fig. 5. Note that for strong interactions, it may be necessary to use realistic interaction 
potentials. It is convenient to express within Eq. (2) by a fitted empirical function 
U2/U = Ai + (A 2 + A 3 a s /a) -\/Vq/E r + A 4 e Asas / a , which describes the on-site interaction 
for two particles but also serves as an approximation for three and four particles [45]. 

8. Mean-field/Gutzwiller phase diagrams 

The SF-MI phase diagram (Fig. 2) for the extended Hamiltonian (1) is derived within mean- 
field theory [2, 4]. By introducing the superfluid order parameter tj) = (6$) = (b\) and by 
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neglecting the fluctuations of quadratic order the lattice sites are decoupled (see Appendix 
F). Within second-order perturbation theory for the Mott state with n particles, we find an 
analytical expression for the SF-MI transition point 

n|J BC (n-l) + l| 2 (n+l)|J BC n+l| 2 

E (n) - E (n - 1) ^ E (n) - E (n + 1) ^ 1 ) 

with J BC = Jnft M °/J^n and the unperturbed energy 

E (n) = {\u™°{n - 1) - v)n/zJ™°. (14) 

Due to the decoupling in the mean-field approach, only tunneling processes with = rij = n 
can be accounted for. Thus, this approach cannot fully implement the occupation-number 
dependent Hamiltonian (1). This can be achieved by using the Gutzwiller approach, where 
the energy functional is minimized with respect to the coefficients f n of the trial wave function 
1^) = Hi J2 n fn \ n )i ( see Appendix G). The Gutzwiller results are depicted in Fig. 2 (red 
lines) and show a nearly perfect agreement with the mean-field calculation (blue lines). The 
mean-field solution (13) can also be applied to the simplified occupation-independent model 
(2) using the lowest-band tunneling parameters J and J B c = Jbc /J as we U as the fitted on- 
site interaction U 2 . Considering the drastic simplifications of this model, the predicted Mott 
lobes (green lines) are in compelling agreement with the results of the occupation-dependent 
Hamiltonian (1). 



9. Conclusions 



We have presented an occupation-dependent model Hamiltonian incorporating multi-orbitally 
renormalized tunneling, density-induced tunneling and on-site interaction. The effectively 
increased tunneling and reduced on-site interaction cause a substantial modification of 
the Bose-Hubbard phase diagram for bosonic atoms. In addition, we have discussed 
a simplified occupation-independent model that includes the lowest-band bond-charge 
tunneling and reproduces the phase diagram well. In general, we have derived a multi-orbital 
renormalization procedure for one- and two-particle processes. It is capable of describing 
strongly-correlated optical lattice systems accurately that cannot be calculated correctly 
by means of variational mean-field treatments [25-28]. Furthermore, it is important to 
include both bond-charge tunneling and multi-orbital corrections as extensions of the Hubbard 
Hamiltonian. Several effects such as the reduction of the bond-charge tunneling cannot 
be described with an effective wave function. The presented results can be transferred to 
fermionic systems, quantum gas mixtures, low-dimensional systems, long-range interactions 
and different lattice geometries. In these systems, the occupation-dependent tunneling may 
also lead to novel quantum phases. Related work on defining a basis with renormalized 
operators can be found in Ref. [46]. 
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Appendix A. Multi-orbital interactions 

The treatment of multi-band Hamiltonians for lattices represents a challenging and still open 
problem posed by its complexity. A common approach is the restriction to the lowest band, 
which neglects interaction-induced admixtures of higher orbitals. In order to account for 
orbital degrees of freedom, we focus on two neighboring sites of the lattice, where we 
apply a multi-orbital diagonalization, so-called configuration interaction method, to each site 
individually (see Fig. Id, main text). In this approach off-site interactions are neglected as 
the occupation of higher orbitals is mainly caused by the local on-site interaction. The three- 
dimensional Wannier functions 

w { L %(r) = w^(r-r L/R ) (A.l) 

constitute the single-particle orbitals, where r L /R = (±§a, 0, 0) is the coordinate of left and 
right site, respectively, a the lattice spacing and 

a = (a xi a y , a z ) (A.2) 

the three-dimensional band index. The maximally localized Wannier functions are obtained 
by means of band- structure calculations for a separable cubic lattice potential V(r) = 
V ^2 • cos 2 (7rr^/a) with the lattice depth V . To include correlations of particles fully, a many- 
particle Fock basis \N) = |n ,ni, ...) with n particles is used, where is the number of 
particles in orbital i. The Hamiltonian for a single lattice site reads 

= J2 e(a)fi{a) + uiafh5) b {a nm^u s \ (a.3) 

where = U a ^U a \ creates and annihilates a particle in the Wannier orbital a 
with single-particle energies e^ a \ In general, the ground state of the single-site problem with 
n particles is given as a superposition of Fock states 

\^(n)) = J2 c N(n)\N(n)) (A.4) 

N 

with real coefficients c N (n). Since the Hamiltonian (A.3) preserves parity the ground state 
consists of many-particle states with even parity. The lowest eigenvalue of the diagonalized 
Hamiltonian matrix directly corresponds to the occupation-dependent on-site interaction 
as discussed in the main text. 
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We model the short-range interaction between the atoms as a box- shaped interaction 
potential 

3 

v{r " r ' } = n °( w - ^ - (A - 5) 

with the Heaviside step function 6 and a variable box width W. The influence of the finite 
interaction range depends only on the ratio Wj a and increases with the lattice depth since the 
wave functions are contracted. However, for W <C a the results (see Fig. 5, main text) are 
almost independent of W. The interaction integrals in Eq. (A.3) are given by 

U^ 5 \g, W) = J d 3 r J Sr 1 w^ a >(r) w^>(r f ) V(r - r') w^\r f ) w^ s \r), (A.6) 

where g is fixed by the condition U^ 000 \g, W) = g J d 3 r \w(r)\ 4 with g = and s- 

wave scattering length a s . Note that the ^-interaction potential can be understood as the 
limiting case of the box interaction potential. However, it behaves differently if including an 
infinite number of orbitals requiring a regularization procedure as described in Ref. [33] (for 
numerical diagonalization see Ref. [47]). For real interaction potentials, the mathematical 
subtlety of the ^-interaction is not present and under certain experimental conditions [35] 
decay processes in higher orbitals may even effectively lead to a finite Hilbert space. Here, for 
the numerical diagonalization nine bands are accounted for in each spacial direction restricted 
to the energetically lowest 6 x 10 3 x n 2 many-particle states. 



Appendix B. Multi-orbital tunneling 

In the following, we consider hopping processes, where one atom tunnels from the left site 
(L) with initially n L particles to the right site (R) with initially n R particles (see Fig. 1(d), 
main text). The tunneling operator accounting for all possible orbital hopping processes reads 
jmo _ J2 af3 Japb^b^ with tunneling matrix elements 

U = -f d 3 r + V(r)) «^(r). (B.l) 

Because of orthogonality relations, only matrix element with a = f3 have non-zero 
contributions. Thus, the tunneling between left and right site simplifies to 

j mo = ^J a ifc )i V£ ) . (b.2) 

a 

with matrix elements J a = J aa , which only depend on the band index a x in x direction. 
The respective tunneling matrix elements J a for the band a x are plotted in Fig. Bl, which 
are positive for even bands and negative for odd bands. The tunneling amplitude is strongly 
enhanced for higher orbitals due to the considerably larger overlap of the Wannier functions. 

As the occupation of higher orbitals is decoupled from off-site interactions and tunneling 
in our model, it is possible to reduce the complex multi-band hopping to an effective hopping 
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Figure B 1 : Single-particle tunneling matrix elements J a for the lowest five bands as a function 
of the lattice depth V and in units of the recoil energy E R = h 2 /8ma 2 . 



J^ nR - It describes the transition from an initial state to the a final state |\P F ) with 

l^i ) = |*K)) L ®|%)> r> 

|*f) = - 1)) L ® |*(n R + 1)) R . 

(see Fig. Id). As discussed in the main text, the effective multi-orbital tunneling can be 
defined as 

JT = -blh JfX (B.4) 

with 

7 mo _ (M i MO |*i) 

™" Vn L (n R + l)' ( 5) 

As inital and final states (B.3) are product states, the expression separates into two terms for 
the left and the right site 

<tf F | J mo |*i> = £ J« (n R + 1)| 6 (Q)t |* (n R )) R 

(B.6) 

x(^(n L -l)|&( Q ) |*(n L )) I 



/L • 



Using the many-particle ground state (A.4) obtained by exact diagonalization, we can define 

j^=(*(n-l)\b^Mn)) 

= CN>(n - 1) c N (n) (N'(n - 1)| S< a > \N(n)) . (R7) 

N'N 



Note that since the coefficients c N can be chosen to be real the conjugated term reads 
fif(n — 1)| U a ^ \*&(n)) = = j^n-i- The effective tunneling can thus be written as 

7 MO _ 1 V 7 7 (q) 7 (q) tB 81 

V n L(n R + 1) « 
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Due to parity conservation of the Hamiltonian (A.3), the ground states at the left and the 
right site are superpositions of Fock states \N(n)) with even parity. Tunneling in an odd 
orbital would alter the parity of both sites resulting in final states with vanishing coefficients. 
Consequently, tunneling can only occur in even orbitals a with positive tunneling matrix 
elements J a . It turns out that the effective tunneling is increased by the tunneling in higher 
bands with the exception of the processes J^ ° = Jf^-i- Here, the lack of interaction in the 
left or the right site prohibits tunneling via higher bands. The interaction at one site, however, 
depopulates the Fock state 0, 0, ...) which results in slightly lowered tunneling energy. The 
process Jf 1 ®, where no interactions are present, is equivalent to the uncorrelated tunneling 
J = </o. 

Intuitively, the mean-field wave functions of the left and right site 

< F (r) = E c n a M Q V) (B.9) 

a 

and 

< F (r) = E c S + i4 Q) W (B.10) 

a 

can be used to calculate the multi-orbital tunneling [28]. The coefficients can be obtained 



from the orbital occupation numbers = y (n^ ) / y/n, with the expectation value 
(nffl) = (*(n) \h^\^(n)). This allows to estimate the multi-orbital by 



« R - / d 3 r < F *(r) (j?- + V(r)j < F (r) 
-Vf( a )f (a) 7 



(B.ll) 



Appendix C. Off-site interactions 

The full lowest-band interaction 

ijkl 

with interaction integrals (3) is commonly restricted to the on-site interaction. However, 
off-site interactions between neighboring sites can reach the same order of magnitude as 
the tunneling. Expanding the lowest-band Hamiltonian accordingly leads to three distinct 
physical processes, namely, bond-charge interaction (density-dependent hopping), correlated 
pair tunneling and density-density interaction [16-19, 37]. While the latter process has to be 
compared with the on-site interaction, the other two constitute tunneling processes. The full 
lowest-band Hamiltonian with nearest-neighbor interaction reads 



(ij) (hi) 

+ \ u ]l ~ 1 ) + v ^ 



2 
3 

(hi) (hi) {hi) 



(hi) 
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Lattice Depth V Q (E R ) 



Figure CI: Lowest-band parameters for on-site interaction U, tunneling J, bond-charge 
tunneling J B c, correlated pair-tunneling J pair and density-density interaction V. 



with matrix elements U = Una, J B c = -Umj 9 </ pa ir = Uujj/2 and V = U^. 

In Fig. CI these parameters are plotted as a function of the lattice depth. The off- 
site density-density interaction is very small compared to the on-site interaction and can 
consequently be neglected. This also applies for the correlated pair-tunneling, which is even 
negligible when compared with the single-particle tunneling matrix element J. The bond- 
charge tunneling matrix element, however, reaches ten percent of the conventional tunneling 
amplitude for intermediate and deep lattices. In addition, it scales with the total particle 
number and can thus be a significant contribution to the total tunneling. Note that the multi- 
orbital renormalization can influence the individual parameters very strongly. 



Appendix D. Multi-orbital bond-charge tunneling 

In analogy to the multi-orbital tunneling, the orbital degree of freedom also affects the bond- 
charge interaction. The single-band bond-charge or density-induced tunneling is introduced 
in the main text [see Eq. (3) and Fig. 1(b)]. As discussed before, we can restrict the calculation 
to two neighboring sites and write the multi-orbital bond-charge operator (see Fig. 1(c), main 
text) as 



tBC,MO _ 1 U a ^f ( T BC _i_ T BC \ £(£)t£W 

J — 2 R V a ^ S J a(3Sj) °L °L 

a(3 7 5 (D.l) 



p ( tbc , tbc n iW)nd) \ m 



with 



tBC 

J a(3~f5 



- Jd 3 r J d 3 r'w£\r)w¥\r') V(r - r') w£\r')w { L S) (r) 
-P J d 3 r J d 3 r' w { "\t)w { 1\t') V{y - r') ^ 7) (r>^(r), 



(D.2) 
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where the Wannier functions are chosen to be real and the sign P = (— i) a *+Ac+7*+<5* i s parity- 
dependent. Note that for a 5-interaction potential, we have = ^/S 7 - The effective 
multi-orbital bond-charge tunneling reads 

tBC,MO ft/- , a 7 BC,MO /t^ q\ 

with occupation-number-dependent matrix elements 

7 BC,MO _ (*f|J BC ' M °|*i) 

V n iA n R + !)( n L + n R - 1) 
Using that and are product states as well as the definitions (B.7) and 

jljftS) = _ i)| S(Wg(7)gW |^( n )) ^ (D.5) 

the effective matrix elements can be computed by 



1 / tBC,MO i tBC,MO\ 

tBC,MO = V 



' .(a) .(/3 7 <5) .(a) .(^) 
Jri L J nR +l Jn R +lJ^L 



^/n L (n R + l)(n L + n R - 1) 



(D.6) 



Note that (*(n - 1)| SWMtSW |$( n )) = £f and that because of the same arguments as 
in the last section, the orbitals a and f3 + 7 + S need to be even; thus the sign P in equation 
(D.l) vanishes. 

Appendix E. Effective potential analogy 

Considering bond-charge tunneling between neighboring sites i and j with the restriction to 
the lowest band, the bond-charge hopping simplifies to 

Jbc = 6j (U m Sft + U w b]b j )b j (E.l) 

(see Eq. (4), main text), where b\ = bf^ and hi = are lowest-band creation and 
annihilation operators, respectively. The interaction integral reads 

U ijk i = g Jd 3 r <(r)^*(r)^(r)^(r) (E.2) 

with lowest-band Wannier functions Wi(r) = w^(r — r^). By inserting the integral 
expressions (E.2) in Eq. (E.l) and using the commutation relations for b\ and 6 i5 the bond- 
charge operator becomes 



Jbc = tybj g (^Jd 3 r w*(r) hi\wi(r)\ 2 Wj(r) 

+ Jd 3 r <(r) (%-l)|^(r)| 2 ^(r) 



(E.3) 
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In order to combine bond-charge hopping and conventional tunneling J = —Jb\bj with 

J=-Jd 3 r <(r) + V(r)^j Wj (r), (E.4) 

we define a reduced density p^- = f^|i^(r)| 2 + (hj — l)\wj(r)\ 2 on neighboring sites i and j 
excluding the hopping particle. Omitting the hopping particle corresponds to the exclusion of 
self-interactions. As a consequence, the total tunneling can be written as 

^ff° = J + Jbc 

(E.5) 



= J <Pr <(r) + V(r) + gp^wfc) b\b r 



In this integral, p^ can be replaced by —> p(r) — \wj\ 2 with the density p(r) = 
Si n d^( r ) | 2 °f a U particles in the lattice, since the Wannier functions are localized. Finally, 
we can identify the effect of the bond-charge interaction as a normal tunneling within the 
effective potential = V(r) + g(p(r) — \wj\ 2 ) experienced by the hopping particle (see 
Fig. 4c, main text). However, since Umj = Uijjj, one may also define p^ = (hi — \)\wi\ 2 + 
(hj — ^)\wj\ 2 . In this case, the effective potential for a homogeneous filling n > 1 can be 
written as V(r) + g^^n — |)|i^| 2 . This form of the potential can be used to perform a 
band structure calculation in the effective potential [21]. Note, however, that in general the 
effective potential is not separable into its spatial directions. Furthermore, in a band structure 
calculus the Wannier functions are adapted to the effective potential whereas in Eq. (E.5) the 
Wannier functions of the bare lattice potential are used. 

Appendix F. Mean-field approach 

The full Hamiltonian including the bond-charge interaction as well as the multi-orbital 
corrections to the tunneling, bond-charge and on-site interactions reads 



!,MO 

(hj) (hj) 



(F.l) 



n 7 , 



(see Eq. 1, main text), where (i, j) denotes nearest-neighbor tunneling and p is the chemical 
potential. The hopping in the renormalized multi-orbital basis is described by b\bj, while hi 
corresponds to the number of particles on site i. In order to apply mean-field theory (see 
Ref. [2, 4]), we introduce the superfluid order parameter ip = (pi) = (b\), where ip ^ 
corresponds to the superfluid phase (SF) and ^ = defines the Mott insulator (MI) with 
a fixed number of atoms per lattice site. The decoupling of the lattice sites is achieved by 
neglecting the fluctuations between b\ and bj of quadratic order, i.e., 

Wbj « Wbj ~ (&! ~ (bjm - fa)) = i>(b\ + 6,-) - ^. (F.2) 
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However, the occupation-number dependent tunneling parameter still couples the sites and it 
is necessary to approximate 

X>® + h) - ^ 2 K% « z I>® + - ^ J £%i> (R3) 

which corresponds to substituting J£±i n . and J%* n .±i by This is necessary, as the 

fluctuations are described by a single parameter i/j 9 which has no internal structure that would 
allow for the implementation of an occupation-dependent tunneling amplitude. This can be 
circumvented by using the Gutzwiller approach, where the information on the population 
of the individual Fock states is preserved (see Appendix G). However, we find that this 
approximation causes only small deviations. Analogously, we neglect quadratic fluctuations 
in the bond-charge term between neighboring lattices sites. Disregarding terms of the order 
C>(^ 3 ) we find 

b\hibj + b\hjbj « ^{b\hi + fijbj). (F.4) 

We perform second-order perturbation theory in of a Mott lobe with n particles per site 
and thus restrict the tunneling to the symmetric terms and «/^' M °. This results in a 
decoupled single- site Hamiltonian Hf s /zJ^° = H + vpV with 

A) = ^Unn(n -l)-ph + 
V = _ft + h)_ftfi + nb)J^ 

z = 6 the number of nearest neighbors, and (U ni p) = (C/M°, J*°' MO , fj)/zJ™°. The 
unperturbed energy for a Mott state \n) = ~^b^ n |0) with n particles is given by 

E Q (n) = (n\ H \n) = ^U n n(n - 1) - fin. (F.6) 
The perturbation series up to third order in %jj reads 

E(V) = E (n) + E 2 {n)^ 2 + 0(^ A ) (F.7) 

with the second-order correction 

F( y n|jg°(n-l) + l| 2 (n + l)|jg°n + l| 2 
2[U) E (n) - E (n - 1) + E (n) - E (n + 1) + ' 

The boundary of the SF-MI phase transition is given by the Landau criterion E 2 (n) = for 
second-order phase transitions, which can be solved for \i (see Fig. 1, main text). 

Appendix G. Gutzwiller approach 

In the mean-field approach elaborated above, only the tunneling matrix elements and 
Jn n ,M ° are taken into account due to the decoupling approximation. Using a Gutzwiller 
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approach [3, 48], the full occupation-dependent Hamiltonian (F.l) can be evaluated without 
this restriction. In this approach, the trial wave function 

oo 

\Gi)=Y,f»\ n )i 

71=0 

with coefficients f n is used, which assumes that the wave function at each site i can be written 
as a sum of local Fock states. In a homogeneous lattice with M sites, all sites are equivalent 
and we can write \G) = rL=i 1^)- The Hamiltonian (F.l) can be split into four parts, namely, 

H = J eff + J^ + U + fi. (G.2) 

The expectation values for this Hamiltonian, depending exclusively on the expectation values 
for two adjacent sites i and j, are given by 

(J eff ) =-zM {G t G 3 \ b^jfX \GiGj) 

= — zM ^T^ fm+lfnifrij+lfnj (G.3) 



(J e B ff c ) = — zM (G^blim + h^X \GiGi) 

= — zM ^ fm+lfnifrij+lfnj (G.4) 
rii ,rij 

x (m + nj)y/ni + ly/nj + lJJ+JJS, 

(U) = l -MY,U n n{n - l)|/ n | 2 , (G.5) 

n 

<A) = -Mj]^n|/ n | 2 , (G.6) 

with z = 6 nearest neighbors. The system is in a Mott-insulator state, when all fluctuations 
vanishes, i.e., when \G) = Yli \ n )i i s ^e energetically lowest state. The resulting Gutzwiller 
phase diagram is plotted together with the mean-field results in Fig. 1 in the main text. Both 
approaches, Gutzwiller and mean-field, agree almost perfectly. 
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